
;; Radius in meters
r = (10.^(findgen(201)/100.+1.)) * 3.08567758d16


;; Mass options
m = (10.^[2.,3.,4.,5.])*1.99d30


;; Temperature options
t = [10.,20.,30.,40.]

sig2 = 3 * 1.38d-23 * t / 2.8 / 1.67d-27

rho = 4.68d-19

alpha = 0.7323d

colors = ['Dodger Blue','Opposite','TG3','Red']
linest = [1,0,3,4]


myps,'./emaf_paper/plots/confinement_pressure.eps',xsize=9,/cmyk
multiplot,[2,1]

plot,[0],[0],/ylog,/xlog,yr=[1d-14,1d-11],/nodata,tit='Jeans Criterion',$
     xtit='R [pc]',ytit='P!dext!n [dyne cm!u-2!n]',xr=[10,1000]

FOR i=0,3 DO BEGIN
   FOR j=0,3 DO BEGIN
      
      pext = rho * (sig2[i] / 3.d - alpha * 6.67d-11 * m[j] / 3. / r) * 10.
      
      cgOplot,r / 3.08567758d16,pext,color=colors[i],linestyle=linest[j],$
              thick=3
      
   ENDFOR
ENDFOR

al_legend,/top,/left,colors=colors,['T = '+string(t,format="(I0)")+' K'],$
          linestyle=0,thick=3,box=0,charsize=0.85
al_legend,/top,/right,box=0,linestyle=linest,thick=3,$
          ['M = 10!u'+string(alog10(m/1.99d30),format="(I0)")+$
           '!n M!d'+cgSymbol('sun')+'!n'],charsize=0.85


multiplot

p0 = 1.6d-12                    ; dyne/cm^2

z = findgen(250)+1              ; pc 

h = 141.                        ; pc

pgal = p0*exp(-z*z/h/h)

plot,z,pgal,xtit='Z [pc]',/ylog,yr=[1d-14,1d-11],thick=3,/xst,$
     tit='Vertical Pressure Profile'


al_legend,box=0,/top,/right,linestyle=0,thick=3,$
          ['P!dext!n = P!d0!n exp[-(z/h)!u2!n]']


myps,/done,/mp


END
